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Abstract 

Grid-based discretizations of the time dependent Schrodinger equation coupled to an external 
magnetic field are converted to manifest gauge invariant discretizations. This is done using general- 
izations of ideas used in classical lattice gauge theory, and the process defined is applicable to a large 
class of discretized differential operators. In particular, popular discretizations such as pseudospec- 
tral discretizations using the fast Fourier transform can be transformed to gauge invariant schemes. 
Also generic gauge invariant versions of generic time integration methods are considered, enabling 
completely gauge invariant calculations of the time dependent Schrodinger equation. Numerical 
examples illuminating the differences between a gauge invariant discretization and conventional 
discretization procedures are also presented. 
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I. INTRODUCTION 



The fundamental laws of physics can (without exceptions) be related to certain continuous 
symmetries. In other words, by requiring that a model should be invariant with respect to 
a certain symmetry, the model is more or less completely determined. The Standard model 
of particle physics and Gravitation |4( are examples of such theories. 

As an example, a model with a complex scalar field, i.e., a model of charged bosons, 
and a requirement of local t/(l)-invariance, or gauge invariance, will immediately yield the 
Maxwell-Klein-Gordon (MKG) theory, which in the non-relativistic limit reduces to Maxwell- 
Schrodinger theory. In addition to defining the theory, the continuous symmetries give rise 
to conserved quantities through Noether's theorem(s) 

aaa 

, and the local [/(l)-symmetry 
of the MKG-model ensures the conservation of local electric charge. 

In particle physics, and especially in the QCD-part of the standard model, numerical 
calculations are often done using Lattice Gauge Theory (LGT) 0,1^, lo|. This is a nu- 
merical procedure, actually motivated from the continuous theory, designed to preserve the 
underlying continuous gauge symmetry. In a previous article this discretization scheme was 
applied to the MKG-equation, with emphasis on the continuous Z7(l)-symmetry and conser- 
vation laws deduced from discrete versions of Noether's theorem(s) ll|. By preserving the 
L^l) -symmetry of the MKG-model on the discrete level, a discrete equivalent of the con- 
servation of local electric charge is immediate, which not only makes the scheme consistent, 
but is also a good indicator of stability. By a more standard discretization of the model, the 
local £7(l)-symmetry is broken, which again implies that the scheme is not consistent with 
the continuous formulation. This will also reveal itself through the fact that the physical 
observables calculated are dependent on the gauge chosen, obviously in conflict with the 
continuous model. 

A similar breaking of the continuous local [/(l)-symmetry has been a known issue with 
discretizations of the Schrodinger equation coupled to an external electromagnetic field. For 
example in atomic physics, results are known to depend on the gauge in which the calculation 
is done - a most unfortunate situation indicating that the calculations are not correct. 
Gauge dependence also leads to interpretation problems of the results. It is the goal of the 
present paper to show how gauge invariant discretizations may be built from existing ones 
with little or no extra effort in the implementations. 
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Simple gauge invariant grid discretizations of Schrodinger operators have been studied in 
previous articles in the LGT formalism with promising results jl3, Q, [l^. The key to suc- 
cess in LGT is that it does not approximate the covariant derivative as a linear combination 
of the gradient and the gauge potential, an element of the Lie-algebra under consideration, 
since such an approximation leads to non-local terms when discretizing the gradient, and a 
question of gauge invariance is meaningless since one compare fields at different spacetime 
points. Instead LGT uses Forward-Euler /central difference approximation of the gradient in 
the various directions, which as argued is not gauge invariant, and then defines the covari- 
ant derivative through the way non-local terms are made gauge invariant in the continuous 
theory. This is done via the Wilson line 18[, to be discussed in the next section, which 
effectively localize non-local terms by parallel transport with the gauge potential as a key in- 
gredient. By defining the covariant derivative in this way, the discrete theory is immediately 
manifestly gauge invariant. 

The aim of this article is to expand the LGT formulation to allow for completely general 
grid discretizations in arbitrary local coordinates of the spatial manifold. Grid discretizations 
are widely used, and include most numerical discretizations of Schrodinger operators in use 
today, such as pseudospectral methods based on the discrete Fourier transform or Chebyshev 
polynomials. We also generalize the discussion to arbitrary coordinate systems, and some 
care is needed in case some of the coordinates are periodic when using global approximations 
(e.g., Chebyshev or Fourier expansions). 

The paper is organized as follows: In Section [III we introduce the time dependent 
Schrodinger equation in general coordinates. In Section II III and [IV] we discuss gauge in- 
variant spatial grid discretizations. We proceed in Section [V] to discuss gauge invariant time 
integration. Finally, in Section IVII we present some numerical results shedding light on the 
difference between gauge invariant and gauge dependent schemes, before we close with some 
concluding remarks in Section IVIII 

II. THE TIME DEPENDENT SCHRODINGER EQUATION AND GAUGE IN- 
VARIANCE 

We consider a particle with charge q and mass m coupled to an external electromagnetic 



(EM) field 16) (E,B). This is a semiclassical approach because the EM- field is obviously 
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affected by the particle, but if we assume that the coupling is weak the approximation can be 
justified. We will work in the non-relativistic regime, but our considerations could easily be 
transmitted to a relativistic model. Moreover, the generalization to more than one particle 
is straightforward, since the EM fields only enter a many-body Hamiltonian at the one-body 
level, i.e., the interparticle interactions are independent of the EM fields. 

We are considering a spacetime domain R x M., with coordinates (t,y), where t 6 I 
is the time coordinate, and y e M. is a point in the spatial domain, usually taken to be 
Euclidean space, but can in general be a Riemannian manifold. In any case, we may work 
in local coordinates x = (x l ), viz, y = y(x) G Ai, with the induced metric tensor gij(x) 
assumed to be time-independent. The wavefunction at some time t is then a complex valued 
scalar function x h- > ifj(x). In addition, the EM-field is described by a gauge potential 
(t,x) I— > (f)(t,x)dt + A(t,x), where is a real valued function and A is a real valued one- 
form. In coordinate basis one usually identifies one-forms with vectors. Thus, if {dx 1 } are 
basis one-forms and {e{\ are basis vectors there is a one-to-one correspondence between 
A = Aidx 1 and A = A l ei. Note, we use the Einstein summation convention except where 
noted. The components of A and A are related by the metric, i.e. A 1 = g^Aj, and the 
physical EM fields (E, B) are given by 

E = -V0 - d t A, B = curl A, (1) 

where we use the shorthand d t = d/dt. 

In the following we work in units such that h — 1. The dynamics of the system is governed 
by the time dependent Schrodinger equation reading 



iD t ijj(t,x) 



~A A + V{t,x) 
2m 



iP(t,x), (2) 



where D t = d t + iqcp is the covariant derivative in the temporal direction, and where the 
"covariant Laplace-Beltrami" operator is defined by 

Aa = —7==Djy^gXxjg jk D k , (3) 

with D k = dk — iqA k being the covariant derivative in the direction k. Moreover, p k = —\D k 
is the (generalized) canonical momentum operator. The term — A a /2m is simply the kinetic 
energy operator, and V(t,x) is an external potential. 
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A fundamental property of the time dependent Schrodinger equation is that it is invariant 
under local gauge transformations, i.e. equation ([2]) is invariant under the following set of 
transformations 



%l)(t,x) i — ^ e iqX{t ' x) ijj(t,x), (4) 
(f)(t,x) i-> (f>(t,x) - d t \{t,x), (5) 
A(t,x) i — ► A(t,x) + VX(t,x), (6) 

where (t, x) \— > X(t, x) is a real valued function meaning that igA(t, x) G u(i), the Lie-algebra 
of U(l) (consult e.g. 0, Q for theory on Lie-groups and Lie-algebras). One says that the 
theory is invariant under local {/(^-transformations, meaning that the physical observables 
are not affected by the transformations. In particular, we note that the electric and magnetic 
fields ([1]) are not affected by the transformations (Q. Moreover, if X A = Xa[{jp^), (x- 7 )] is an 
observable, then the expectation value (if), X A ip) is gauge invariant, viz, 

(e i9 >, X A+vx e^) = (V, X A iP), (7) 

where (•, •) denotes the standard inner product in L 2 (M). 
The usual way to write Eqn. (j2J) is 



id t ijj{t,x) = H{t)^{t,x) (8) 

ij)(t,x), 



-—A A + V(t,x) + q(j)(t,x) 
2m 



where the Hamiltonian H(t) depends on the fields (A, 0). It is well-known that for any two 
t,t', the formal solution to (jHJ) is given by ip{t,x) = U(t,t')il)(t,x f ), where the propagator U 
is 

(9) 



U(t,t')=Texp ^-ijf H(s)ds^j 



with T being the standard time-ordering operator. The propagator depends on the fields 
(A, <j)) in the case of the current Hamiltonian, and under a gauge transformation with pa- 
rameter \(t, x) we have 

U A >A*, = e iqX{t) U A At, t')e- i9A(4,) , (10) 
where A' = A + VA and 6' = 6 - d t X. 
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(r,e)eGc 



{y\f)=y{G)<lJZ 




FIG. 1: A grid in polar coordinates. In local coordinates (r, 0) G M 2 the grid G is a Cartesian 
product, while in on the manifold A4 C M. 2 the grid instead has a (discrete) rotational symmetry. 
A coordinate curve for constant r is also illustrated. In this case, the curve becomes a circle. 

III. DISCRETIZATION ON A SPATIAL GRID 

Many discretizations of Eqn. ([2]) approximate the wave function ip(t,x) at a given time 
t G R on a finite grid G (i.e., y(G) C .M) in order to obtain a semi-discrete formulation in 
which if)(t,-) G L 2 [y(G)] still depends continuously on time. We shall here consider grids 
which in the local coordinates are Cartesian products of one-dimensional grids, i.e., 

G = G 1 x G 2 x ■•• x G n , (11) 

where 

G 3 = { X li X 2i ' ' ' i X Nj }' X k < X k+1- (12) 

Figure [T] illustrates this in the case of polar coordinates in the plane. 
We may list the elements of G using multi-indices, i.e., 

G = {x a = «, ■ • • ,x n a J : Vj, 0<aj< N,} , (13) 

and the multi-indices may again be mapped one-to-one with {0, 1, • • • , N — 1}, where iV = 
N1N2 ■ ■ • N n is the total number of grid points. Thus, we obtain a discrete Hilbert space 
H(G) ~ L 2 [y(G)} of dimension N. 

The natural basis to use in the space 7i{G) is the set of functions e a such that e a (xp) = 
5 a p. These functions are referred to as the cardinal basis [3] or the nodal basis. For any 
ip G TC(G), we now have 

i> = ^2^( x a)e a . (14) 
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A linear operator X onH may be represented by its action on this basis, which determines 
aniVxiV matrix with elements X al 3, viz, 

X aP =[Xe a \{xp). (15) 

Thus, for V G H(G), 

[X^}(x a ) = Y J X^(x p ). (16) 

At times, we will omit the brackets and write Xif)(x a ) for the product Xip evaluated at x a , 
as usually there is no danger of confusion. Likewise, multiplication by u G 7i(G) (or any 
continuous u over M) defines a linear operator, and for ip e H(G) we will write this simply 
as uip. 

The grid discretization invariably comes with a discrete approximation to the (field-free, 
i.e., A = 0) Laplace-Beltrami operator A , although the exact procedure to fix this discrete 
operator may vary. We may assume that this discretization is composed of discrete deriva- 
tives in each spatial direction x j composed with some fixed functions of the coordinates, but 
the exact form is of no consequence to us. 

To clarify these statements, consider for example Cartesian coordinates in two spatial 
dimensions for which A = d 2 + d 2 . In a finite difference approximation we typically have 
a standard 5-point central difference stencil, i.e., 

A ip(x\ x 2 ) « -(<*}<*! + 5lS 2 )ip(x\ x 2 ), (17) 

where Sj is a forward difference, — 5j a backward difference, so that — S^5j is the standard 
3-point central difference operator in the x J -direction, viz, 

- 5}5jf(xi) = ±;[f(xt + h)- 2/(^) + f(x* - h)l (18) 

with h being the mesh width. We have suppressed other spatial coordinates than Xj in the 
latter equation. 

As a different example, consider polar coordinates (a; 1 , a; 2 ) = (r,9) in two dimensions, for 
which 

A = -d r rd r + \d 2 . (19) 

The form of the radial part of this operator leads initially to several different schemes by 
either expressing it as 

-d r rd r = d 2 . + -d r (20) 
r r 
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and then discretizing, or instead attack the original difference operator. In a general coor- 
dinate system there will of course be even more possibilities. 

In any case, the discrete Laplace-Beltrami operator will be on the form 

A «A 0)h = Ao )h (« jl i} ) i) ) (21) 

where 6* are arbitrary approximations to each partial derivative dj. We abuse notation 
a little, as in general we allow 8j ^ ($j) k - m ^he Cartesian coordinate example above, 
5 2 = — 8\8a. In general, however, we assume that like in this example, 5^ is the product of 
k discrete derivative operators Sj/, 1 < £ < k, so that 

A « A 0ih = A 0ih (5 jifc ,x), (22) 

with (5j 5 o = <5j- 

The usual way to discretize A^, on the other hand, which we here will call a "naive" 
discretization, is to employ a similar "recipe" as in the examples to Eqn. ([3]) after insertion 
of Dj = dj — iqAj and simplifying the expression. This, however, always leads to non-gauge 
invariant discretizations, as we will discuss in Section IIVI 

As an example of the naive approach, consider again the polar coordinate case, and for 

simplicity assume A r = for which we obtain 

1 1 
Aa = ~d r rd r + —(d e - iqA e ) 2 

1 a 2 
= A - iq^(d e A e + A e d e ) - ^Aj. 

Assuming further that dgAg = 0, i.e., that A is given in the Coulomb gauge, we get 

A A = A - i2q^A e d e - ^A 2 g . (23) 

The naive discretization of Eqn. fl23|) is then given by inserting the usual grid discretizations 
of d r , d 2 and d 2 . 

Our prescription for a manifestly gauge invariant discretization of A^ in Eqn. ([3]) is simply 
to replace all occurrences of approximations 5j± of dj in the field-free naive discretization 
with a certain corresponding approximation Djk to Dj, derived using methods from LGT 
as mentioned in the Introduction, and whose final expression is given in Eqn. ( 1491) below. 
In other words, 

( 24 ) 
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which will be gauge invariant. This approximation is often quite different from the standard 
naive discretization. 



IV. DEFINITION OF THE DISCRETE CO VARIANT DERIVATIVE 
A. One dimensional manifolds 

1. Gauge transformations 

Consider first the case when M. is a one- dimensional manifold M. = M. 1 d W 1 . The 
reason for this is that in the general case, the differentiation operator dj can be viewed 
as a differential operator on the coordinate curves, these being one-dimensional manifolds. 
Similarly, a generic discrete Sj can be viewed as an operator acting on grid functions over 
the one-dimensional "coordinate grids" obtained by fixing all but the j'th component of the 
multi-index a in Eqn. (fl3l) . Equivalently, Sj defines a discrete differentiation operator acting 
on functions over a discretization of the coordinate curve; see Fig. HJ 

Any one- dimensional manifold M. will either be topologically equivalent to a circle or an 
interval, which may be bounded or unbounded. For example, in polar coordinates (x 1 , x 2 ) = 
(r, 9) in M 2 the angular coordinate curves are circles of radius r while the radial coordinate 
curves are rays from the origin r = to infinity with an angle 9 relative to the x-axis. 

We write x = x 1 for the sole coordinate, omit the time dependence, and Da = d x — iqA(x) 
for the covariant derivative. Under gauge transformations, Da transforms as 



where X'(x) = d x X(x). Intuitively, since M. is one-dimensional, one should be able to 
transform away A(x) completely, by selecting A' = —A. However, if M. is (topologically) a 
circle (with the point x = identified with x = L, for simplicity), this is not possible: There 
are one-forms A(x) which are not the derivative of some zero-form X(x). On the circle, it is 
precisely the constant functions A(x) = A , since then X(x) = A x + b is not a zero-form: 
it is not periodic in x unless Aq = 0! If, on the other hand, M. is topologically an interval, 
A(x) may be transformed away. 

These considerations may become clearer when we observe that, locally, we may write 




(25) 




(26) 
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where 

u(x) = exp ^-iq J A(s) ds^j . (27) 

Whenever M. is topologically an interval, we can choose A' = —A, and this expression is 
global, since then 

u(x) = exp(igA(x)). (28) 
For M. being topologically a circle, no such X(x) exists globally, unless A = 0. 

2. Local approximations 

Let T~L{G) be the discrete Hilbert space corresponding to an iV-point discretization of 
M. = M 1 , being either a circle or an interval as described above. Thus, G C M 1 is given 
by 

G = {x!,x 2 , • ■ -x N }, x k <x k+1 . (29) 

In the case of a circle, we identify x^v+i an d X\ to impose periodic boundary conditions. We 
let h = min(ajfc, iCfe+i) be the mesh width, and typically h ~ 1/N. 

Let 6h be a discrete differential operator on 7i(N), and we assume for the moment that 
5h is a local operator, in the sense that as h — > 0, only a finite number of points in the 
neighborhood of x k € G are used to differentiate ip(xk). As a consequence, there is a largest 
p > such that for any smooth function ip(x) over 

= + (30) 

where the term 0(h p ) is equal to the truncation error, and we say that 5h is a p'th order 
approximation. Examples of local discretizations are finite differences of any order, but not 
pseudospectral methods using for example Chebyshev polynomials or the discrete Fourier 
transform. 

For any (discrete or smooth) ip, the naive discretization t>A,h of Da reads 

D A ,h^(x k ) = (Sh ~ \qA)i)(x k ). (31) 

This operator is not gauge invariant. Let \{x) be given, and consider 

(D A+x ,, h e iX i/>)(x k ) = {8 h e^){x k ) -iq{Ae^){x k ) 
£ e^D A , h ^(x k ). 
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Clearly, this comes about since 

(4e igA )(x fc ) = iq\'(x k ) exp(iq\(x k )) + 0(h p ) 

is only an approximation. 

However, the continuous covariant derivative Da comes about if one tries to construct 

I I 

a gauge invariant classical field theory [I]: since for any differentiate ip, ip{x) and ip(y) 
transforms differently if x ^ y, the limit 

limhifj(x + h)-tfj(x)} (32) 
h— >o h 

has no simple transformation law. Notice that for finite h, [ip(x+h)—ip(x)}/h is the standard 
forward difference operator. We could just as well consider the limit 

\im8hip(x) (33) 

for any local discrete differentiation operator. 

Introducing a comparator function U(x, y) with the transformation law 

U(x, y) — ► e iqX(x) U(x, y)e" i<?Afe) , (34) 

we see that for any y, the function U(x, y)ip{y) transforms in the same way as ip(x). Explic- 
itly, the comparator is given by 

U(x,y) = e~ Ul % A ® dt . (35) 

For any finite h, consider again the discrete difference operator 5h applied to U(x,y)^(y), 
but acting on the variable y, i.e., 

D h ip(x k ) = [8 hi yU(x k ,y)ij;(y)](x k ). (36) 

The notation implies that the discrete derivative is evaluated ai y = x k - This operator is 
obviously gauge invariant, so that 

Da+M^WM = e iX D A ^(x k ). (37) 

The path from x to y is in general ambiguous if M. is a circle: we may move either 
clockwise or anti-clockwise, and also through several revolutions before ending up at y. 
However, as h — > 0, we desire our discrete Da,k to converge to Da- The truncation error 
does not vanish unless we choose the shortest path, with vanishing length. Then, Eqn. (130]) 
implies that for any smooth i/j(x), 

[D A ,hip](x) = D A ^{x) + 0(h p ). (38) 
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3. Global approximations 



The fact that 8h was a local approximation to the derivative was crucial above, as it 
allowed us to resolve path ambiguity. For a global approximation this is not the case. 

A global approximation to d x in general has exponential order of approximation as it 
utilizes all grid points Xk G G to estimate the derivative. That is, for any smooth i/j(x), the 
truncation error is 0(h N ) = 0(/i 1///l ), i.e. 

d h ip(x) = d x ^(x) + 0(h N ), (39) 

so that the order of approximation in fact increases as h —>■ 0. 

As U(x, y) may depend on path, and as x and y have arbitrary separation in a global 
method, the limit h — > does not resolve the path ambiguity. 

The only way to overcome this, is to ensure that the comparator itself is path- independent. 
This is the case if and only if A(x) = d x A(x) for some smooth A, since then the fundamental 
theorem of analysis yields 

f A(t)dt = A(x)-A(x'), (40) 

Jx' 

independently of the path taken. For M. being a circle, this means that A(x) must be the 
derivative of a periodic function. 

We therefore decompose the one-form A(x) as 

A(x) = A + A 1 (x), (41) 

where A is a constant such that Ai(x) = d x A(x) for some smooth X(x). Formally, A is 
the projection of A(x) onto the orthogonal complement of the range d x , i.e., Ran(9 a; )- L . The 
decomposition (T4~T1) is unique and it always exists. We may say that A\ is the "largest part 
of A that may be transformed away." We then obtain 

D A = e i ^ XA ^d x -iqA )e- i ^ XA ^ (42) 

for the covariant derivative. It is clear that if Aq is nonzero, it may never be transformed 
away using a gauge transformation. 

In the case of M. being an interval, Ran(9 x ) ± = {0} , implying A Q = 0, and we get 

D A = e nS*A{t)dt dxe -i q S*A{t)dt^ (43) 
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where any anti-derivative of A(x) may be used. For Ai being a circle, however, 

Ran (d x ) ± = { constant functions }, (44) 
since J Aq = Aqx + b is not periodic unless Aq = 0. It is straightforward to show that 

A = (A) = y [ L A{t)dt. (45) 
L Jo 



We define a modified path independent comparator given by 



U (x, y) = exp 



iq / (A(s) - Aq) ds 



(46) 



Since it is independent of path, we may write 

U(x,y) = u*(x)u(y), u(x) = U(x ,x), (47) 

where xq G G is any reference point. Combining Eqns. f|42l) and (|36|) we get 

£>A,hip(xk) = [bh, y U (x k ,y)ip(y)](x k ) - iqA tp(x k ), (48) 

and using Eqn. (1471) we may rewrite this as 

D A ,h^{xk) = u*5 h uip{x k ) - iqA tp(x k ), (49) 

which may be a more practical expression to implement. This covariant derivative is valid 
for any one-dimensional manifold topology and any discretization of the derivative, and we 
note that in particular for A = 0, it is equivalent to the original expression (1361) for local 
discrete derivatives. 

B. General manifolds 

For global methods, the one-dimensional case necessitated the computation of Aq given 
by Eqn. fj45|) . As the covariant derivative in this case was gauge equivalent to using the naive 
discretization with a constant A(x) = A , one may wonder what we have to gain from the 
approach in this case: Why not use the standard naive discretization using this particular, 
and physically equivalent, gauge? Most manifolds are, however, not one- dimensional. In 
this section, the case of M. = M. n being a general n-dimensional manifold is treated by 
simply defining Da,k for each spatial direction. In this case it is in general not true that 
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the method becomes gauge equivalent to a naive discretization: we may not find an A such 
that the problem may be solved with a naive discretization. Said in another way, on Jvt n , 
the splitting (j4"Tl) becomes 

A(x) = A (x) +VA(x), (50) 

where Aq G Ran(V)" 1 is the part of A which may not be transformed away, and this is of 
course not a constant function in general. 

In the z'th direction, at the point y(x) G Jvt, the continuous covariant derivative is given 

by 

D(A)i = d i -iqA i (x), (51) 

being an operator that constructs the z'th component of a one-form field, i.e., D(A)if)(x) = 
[V — iqA(x)]ip(x) is a one-form field with components D(A)iijj(x). 

As discussed in Section IHIt we are given discretized derivative operators 5i (we suppress 
the subscript "/t" in the sequel, and also the distinction between local and global discrete 
differentiation operators di) which only involves grid points along the z'th coordinate curve 
at x a G G = G 1 x G 2 x ■ ■ • x G n ; see Fig. [TJ for an illustration. Thus, <5j may be viewed 
as a discrete derivative on discretization of a one- dimensional manifold (the z'th coordinate 
curve at x a ) with grid G l . From Section HVAl we then have the discrete covariant derivatives 
D(A)i given by 

D(A)i = u*(x)5 iUi (x) - iqA (x) h (52) 

where A (x)i is the quantity A in Eqn. (j4"Tl) for D(A)i - in general not a constant since 
it depends on the other coordinates x\ j ^ i. Neither is it given by the decomposition 
( |50l) . Moreover, Ui(x) (or more precisely u*(x)iii(y)) is the corresponding path- independent 
comparator for differentiation in the z'th direction. 

To be precise, we write out these quantities in the general case. 

The quantity A (x a )i is zero for coordinate curves that are topological intervals, such 
as the radial coordinate curves in polar coordinates. For periodic coordinates, such as the 
angle 9 in polar coordinates, the coordinate curves are circles. In that case, 

1 f Li 

A (x a ) l = -—— Ai(x\--- ,x t ~ 1 ,s,x l+1 ,---) ds, (53) 

J^i{X a ) JO 

where Lj is the length of the coordinate curve. In polar coordinates, x a = (r a ,9 a ), and 
Li = 2nr a . 
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The comparator function Ui(x a ) is now given by 



Ui(x a ) = exp 



-iq (Ai(- ■ ■ ,s, ■ • ■) - A (x a )) ds 




(54) 



where the arbitrary reference point has been chosen as x l a = 0. 

Gauge invariance of D(A)i follows from the gauge invariance in the one-dimensional 
case. Clearly, gauge invariance necessitates calculating the comparator functions and Ao(x)i. 
However, these enter the discretizations only as multiplicative operators which are diagonal 
in the nodal basis. It is therefore a one-time calculation inducing little overhead in general. 

V. TIME DISCRETIZATION 

Thus far we have studied a discretization of the time-dependent Schrodinger equation in 
continuous time. However, when solving the problem numerically one needs to discretize 
the model in time as well. Again, with inspiration from classical LGT this can be done 
manifestly gauge invariant for every scheme with a grid based approximation of the time 
derivative. 

The standard way to propagate the Schrodinger equation (j2J) is to attack the form (jSJ) 
instead, viz, 



and then use standard techniques to integrate, analogously to the naive spatial discretiza- 
tions. However, this will of course lead to non-gauge invariant solutions. 

Let us consider how gauge invariant formulations of some simple schemes can be con- 
structed. For simplicity, we will assume that the wave function ip(t,x) is only sampled at 
equally spaced points in time, i.e., t n = nr with n = 0, 1, . . .. At each time t n , we write 
ip n G 7i{N) for the corresponding spatially discrete wave function. 

Let 5 t = dt + 0(r k ) be a local approximation, and assume that a naive discretization of 
a generic Schrodinger equation with Hamiltonian H(t) is given by 



\j\<m 

where H n is the Hamiltonian H(t) evaluated at t — nr. Here, Cj are constants, the notation 
indicating that only a finite number of such constants are involved. Such schemes include 




(55) 




(56) 
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the standard implicit Crank-Nicholson and leap-frog schemes [19j. For the Crank-Nicholson 
scheme, St is the forward Euler discretization, while cq — c\ — 1/2 (m = 1 and c_i = 0). 
For the Leap-Frog scheme, St is the centered difference with step length 2r and Cq = 1 (with 
m = 0). Using similar considerations as in Section HVl for spatial differentiation operators, 
the corresponding gauge invariant discretization of d2J) for arbitrary fields becomes 

ib t r = £ t n ,H n A +j r +j , (57) 

\j\<m 

where c n j = CjU(t n ,t n+ j), with 

17(<, *') = exp (i jf" 0(s)rfsj (58) 

being the comparator for the time coordinate. The Hamiltonian H^(t) is given by 

#AW = -7^A A + U(t,x), (59) 

and excludes the term qcj)(t,x) which is now absorbed into the covariant derivative D t . 

We now observe something peculiar: The scalar field <p(t, x) may be transformed away 
globally by the gauge parameter A = f (f)(s)ds, yielding the so-called temporal gauge. In 
this gauge U(t,t') = 1, so the gauge invariant scheme fl57|) reduces to the naive scheme f )56|) 
- of course with a different Hamiltonian Ha+vx- 

In fact, these considerations hold for any gauge invariant numerical integration scheme: 
gauge invariance implies that the temporal gauge in particular may be used, for which the 
integration method reduces to the naive non-gauge invariant scheme applied to Ha+vx- 
Notice, however, that the latter operator is time dependent, even if V and the fields A and 
(p are time-independent functions. 

To make this statement precise, let Uh(t + r, t) be a general numerical propagation scheme 
for a generic Hamiltonian H(t), i.e., it approximates the propagator U(t + r,t) in Eqn. 0, 
viz, 

U(t + r,t) =U h (t + r,t) + 0(r m ), (60) 

where 0(r m ) is the truncation error of the scheme. Thus, the wave function if) n is propagated 
by 

^P n+1 =U h (t + r,t)^ n . (61) 
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Assuming that Ua,4> is a gauge invariant generalization of Uh applied to the Hamiltonian 
Ha, it must transform according to Eqn. (flOl) . The temporal gauge is achieved by selecting 
A = A given by 

A(t,x) = / 4>(x,s) ds, (62) 
Jo 

which gives 

A' = A + X7A(t,x) = A + X7 [ <f>(x,s)ds. (63) 

Jo 

We obtain 

U A<t> (t, = e i?A ^ +V A,o(t, t')e-^\ (64) 

where Ua+va,o must be equal to the original gauge dependent propagator applied to the 
Hamiltonian Ha+va- 

It is not always easy to identify an expression for UA,<p{t, t') in a general gauge, but from 
the above considerations, the temporal gauge is sufficient anyway. 

The selection of a particular gauge for time integration may seem unnatural. However, it 
is the structure of the Schrodinger equation together with the fact that any field (j) may be 
transformed away that yields this conclusion. The vector potential A cannot in general be 
transformed away - therefore we should not choose a particular gauge for spatial operators. 



VI. NUMERICAL EXAMPLE 



In Ref. [13J, some promising gauge invariant eigenvalue calculations are shown using the 



classical LGT formalism, i.e., with standard finite differences in space. In Ref. [15] higher- 
order finite differences are used, and the results are equally promising. Even though the 
published experiments are all with the "standard" example of a uniform magnetic field in 
the z-direction applied to a planar system, there is little doubt that the gauge invariant 
formulations offer favorable properties over the non-gauge invariant methods, as there are 
always gauges that behave very badly. One may simply choose a rapidly oscillating gauge 
parameter X(t,x) to completely destroy the accuracy. Choosing the "right gauge" in a non- 
gauge invariant scheme may not at all be simple or even possible. In any case, a gauge 
independent method will "factor out" any non-physical effect of the choice of gauge, making 
interpretations easier, and it is reasonable to expect that gauge invariance should stabilize 
the discretization because of this. 
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We focus here on time integration only. The benefit of employing spatially gauge invariant 
schemes has already been established in e.g. 13)]. A gauge invariant discretization of the 
time-dependent Schrodinger equation could enable practitioners to push the limits of what 
is possible to compute and interpret. 

We consider a single particle in a one-dimensional system; a very simple system but one 
whose numerical properties are reflected in more realistic settings. We set m = q = 1, and 
consider the Schrodinger equation ([2]) on the form 

i[d t + i4>(t,x)]i;(t,x) = ~[d x -iA(t,x)] 2 i/;(t,x). (65) 

We consider a spatial truncation [-L/2, +L/2] C M, and use a finite difference discretization 
with N + 2 equally spaced grid points x k = kh E G, k = 0, 1, . . . , iV + 1. The grid spacing is 
given by h = L/(N + 1). Thus, at a time t, if)(t,x) ~ ip(t,Xj) G 'HiG) is our discrete wave 
function. 

A common situation in atomic physics arise when one considers the so-called dipole 
approximation [jjj], in which the fields A and (j) take the form 

<P(t,x) = f(t)x 

A(t,x) = 0, (66) 

corresponding to a time-dependent electric field E(i) = —f(t)e x and B = 0. These fields 
are of course not solutions of Maxwell's equations. The particular gauge in Eqn. ( 1661) is 
referred to as the length-gauge. The so-called velocity gauge is obtained by transforming 
away 4>(t,x) using the gauge parameter X(t,x) = f* 0(s, x)ds, i.e., it is the temporal gauge. 
We obtain 

<i>\t,x) = 0, 

A'{t,x) = [ d x (p(s,x)ds. (67) 
Jo 

The wave functions in the two gauges are of course related by ip'(t, x) = e lX ^' x ^ip(t, x). These 
two gauges are commonly studied, and may give different results in actual calculations; a 
sure sign of a significant error. 
A common choice for f(t) is 

f(t) = csm(ut), (68) 
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describing an oscillating electric field with frequency lo. 

Using finite differences, a typical naive semi-discretization of Eqn. fl65|) is 



Sdrf&xj) = -~(6 + -L4)(<T-L4) +4> ^(t,Xj) 

= ~ [-5+5- + i(A5~ + 5 + A) + A 2 + 20] j>(t, Xj ). 



(69) 



where 5 + is a forward difference, and 5" is a backward difference. As earlier, A and should 
be interpreted as diagonal multiplication operators. As = — S~, it is easy to see that 

the operator on the right hand side of Eqn. fl69l) is actually Hermitian. 

The corresponding gauge invariant semi-discretization is, in the temporal gauge, 



with u(t,x) = exp[— i\(t, x) + iA(i, 0)] being the comparator function. 

To integrate Eqns. ( 1691) and ( 1701) in time, we select a somewhat non-standard approach. 
It is well-known that an approximation to Uit + r,t) for a given Hamiltonian H{t) is 



where the error is 0{j 2 ). Propagating ip n = if)(t n ,Xj) G T~t(G) using U T {t) gives an error 
increasing roughly linearly as function of the number of time steps. The integral is evaluated 
using Gauss-Legendre quadrature using two evaluation points, giving practically no error in 
the integral as long as fit) does not oscillate too rapidly. 

We choose c = 5 and u = 10 for the electrical field, and integrate for t n = nr < 6tc/u, so 
that the electric field oscillate exactly three times before we terminate the calculation. The 
spatial domain is of length L = 40, and we use iV = 255 points. 

We choose ip(x, 0) = exp(— x 2 /2) as initial condition (which is normalized numerically in 
the calculations). The analytic solution using this particular problem (on whole of K) can 
be computed in closed form, but we choose instead to perform a reference calculation using 
a pseudospectral discretization using N + 1 points and a much smaller time step, giving in 
this case practically no error. 

Figure [2] shows the error \\ip n — ^ ex act(^n)|| as function of t in the three cases. Clearly, the 
velocity gauge has somewhat smaller error, and also the length gauge and gauge invariant 
calculations have almost indistinguishable errors. The latter fact can easily be understood 




(70) 



U{t + T,t)n Ur(t) = e-*Ji TH(s)ds , 



(71) 
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FIG. 2: (Color online) Errors in the time-integration for the electric field E(t,x) = — c$in(ut)e x . 
The velocity gauge (dashed/dotted green) has somewhat larger error than the gauge invariant 
(solid black) and the length gauge (dotted red) calculations. The low error in the velocity gauge 
is a "stroke of luck" when compared with Fig. O where the velocity gauge has the largest error. 
Notice that the oscillations of the EM-fields clearly affect the errors. 

by inserting ip' = exp(i\)ip into the semi-discrete formulations, and noting that 5 + 5~ and 
u*5 + 5~u are unitarily equivalent, i.e., having the same eigenvalues. The semidiscrete equa- 
tions are thus actually equivalent, and any discrepancy showing in the graphs for the gauge 
invariant scheme and the length-gauge calculation comes from errors in the time integration. 

The spectrum of (d~ + — iA)(5~ — iA) is not equivalent to that of 5 + 5~ when A / 0, 
however. In fact, it is readily established that the latter operator has eigenvalues depending 
strongly on A, and therefore on the particular gauge used. Hence, it is expected that the 
velocity gauge, or any other gauge in which A ^ 0, should perform worse than either of the 
other gauges in the generic case. 

To test this statement, we perform a calculation using a different field <p(t, x) in the length 
gauge, namely 

4>(t,x) = csm(u;t — fix), (72) 

which may describe incoming electromagnetic waves (e.g., a laser) along the x axis. Figure 
|3] shows the errors as function of t in this case, using fi = 1, clearly showing that the 
velocity gauge indeed has the larger error. Moreover, a "mixed gauge" calculation is shown, 
where the length gauge fields are transformed using a gauge parameter X(t,x) = x, chosen 
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FIG. 3: (Color online) Errors in the time-integration for the electric field E(t,x) = c/xcos(wi — 
fix)e x . The gauge invariant (solid black) and the length gauge (dotted red) calculations have the 
smallest errors, while the velocity gauge (dot-dash, green) has clearly the largest error. This should 
be contrasted with the results in Fig. [21 where the velocity gauge has the smallest error. 

somewhat arbitrarily. Now, both A and are non- vanishing, and the error is seen to behave 
accordingly. 

We notice that the gauge invariant calculations in both cases are well-behaved, and no 
choice of gauge will of course affect the calculations. Moreover, the length gauge is equivalent 
to the gauge invariant calculations only when A = 0; this holds in general in one dimensional 
systems, but of course not in arbitrary dimensions, where the magnetic field usually cannot 
be transformed away in this way. 

VII. CONCLUSION 

We have discussed a method based on LGT to convert virtually any grid-based scheme 
for the time- dependent Schrodinger equation (jSj) to a gauge invariant scheme, in both space 
and time. We have considered discretization in arbitrary coordinates on arbitrary spatial 
manifolds. The theory is directly generalizable to many-particle systems as the EM-fields 
obviously only enter at one-body level in the many-body Hamiltonian. Moreover, the compu- 
tational overhead of the gauge invariant schemes compared to the original ones are negligible. 

Our numerical simulations of time-dependent problems, albeit simplistic, indicate that 
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the gauge invariant schemes perform on average better than standard schemes, even though 
the original "naive" scheme may be better in specific cases. 

A further line of work would be to rigorously understand the accuracy gained by intro- 
ducing gauge invariance. 
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